MIT-CTP 3964 



Solitonic ground states in (color-) superconductivity 

Dominik Nickel 1 and Michael Buballa 2 

'Center for Theoretical Physics, MIT, Cambridge, MA 02139, USA 
2 Institut fur Kernphysik, Technische Universitat Darmstadt, Germany 
(Dated: March 5, 2009) 

We present a general framework for analyzing inhomogeneous (color-) superconducting phases 
in mean-field approximation without restriction to the Ginzburg-Landau approach. As a first ap- 
plication, we calculate real gap functions with general one-dimensional periodic structures for a 
3 + 1-dimensional toy model having two fermion species. The resulting solutions are energetically 
favored against homogeneous superconducting (BCS) and normal conducting phases in a window for 
the chemical potential difference Sfi which is about twice as wide as for the most simple plane- wave 
ansatz ( "Fulde-Ferrell phase" ) . At the lower end of this window, we observe the formation of a soli- 
ton lattice and a continuous phase transition to the BCS phase. At the higher end of the window the 
gap functions are sinusoidal, and the transition to the normal conducting phase is of first order. We 
also discuss the quasiparticle excitation spectrum in the inhomogeneous phase. Finally, we compare 
the gap functions with the known analytical solutions of the 1 + 1-dimensional theory. 



I. INTRODUCTION 

Inhomogeneous ground states due to imbalanced Fermi 
surfaces have been discussed in various contexts. Theo- 
retical investigations started off by considering a clean 
paramagnetic superconductor exposed to an external 
magnetic field. For such a system Fulde and Ferrell an- 
alyzed the ground state with the order parameter, i.e., 
the gap function, forming a plane wave Larkin and 
Ovchinnikov extended their work by considering more 
general inhomogeneous ground states, but relying on the 
Ginzburg-Landau expansion 

In recent years these ideas have attracted new inter- 
est from two fields. One of them are systems of ultra- 
cold atoms [3(, where new fascinating techniques open 
unprecedented possibilities to study the pairing of im- 
balanced Fermi systems in a trap. Here, unlike in solids 
where the electron interaction is often difficult do under- 
stand in detail, imbalanced systems can be prepared in a 
rather straight forward and controlled way. 

In this paper, we will mainly aim at color supercon- 
ductivity in deconfined quark matter. Here the problem 
of imbalanced Fermi surfaces is almost unavoidable. It is 
expected that color superconducting phases are present 
in the QCD phase diagram at sufficiently high densities 
and low temperatures (see Refs. 0, H, S 0, II, H 03 for 
corresponding reviews). In nature, the most promising 
places to find these conditions are the centers of neu- 
tron stars. Here the system must be in beta equilibrium 
and, at least globally, electrically and color neutral. This 
would be fulfilled if there were equally many up-, down-, 
and strange quarks. However, at densities which can be 
reached in neutron stars, strange quarks are expected to 
be considerably suppressed by their mass. In turn, this 
forces the density of down quarks to be larger than the 
density of up quarks in order to achieve electric neutral- 
ity. Since, on the other hand, the most attractive channels 



involve quarks of unequal flavors, we are naturally led to 
the problem of pairin g in an imbalanced Fermi system 
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Let us briefly recall what the problem actually is. In 
BCS theory, pairing occurs among fermions with oppo- 
site momenta, forming Cooper pairs with zero total mo- 
mentum. If both fermions are at their respective Fermi 
surface, the pair can be created at no free-energy cost and 
the pairing is is always favored as soon as there is an at- 
tractive interaction. This is, however, no longer the case 
if the Fermi momenta of the fermions to be paired are un- 
equal. BCS pairing then requires that the Fermi spheres 
first have to be equalized. In the case of quark matter this 
could be realized, e.g., in a weak process which replaces 
some of the down quarks by strange quarks. Of course, 
this will only be favorable if the free energy which is 
needed for this process is overcompensated by the pair- 
ing energy. This sets a limit for this mechanism in terms 
of the Fermi momentum difference in the unpaired sys- 
tem and the BCS gap 

Therefore the question arises how the system reacts 
if the imbalance does no longer allow for BCS-like pair- 
ing. Sticking to homogeneous phases, some authors have 
suggest ed so-called gapless or breached pairing phases 
[23, |2J, |25|, |26| , where equal Fermi surfaces are created by 
lifting some of the fermions to higher momentum states. 
At these "new" Fermi surfaces the fermions can again 
form Cooper pairs with zero total momentum. It was 
found, however, that this pairing mechanism suffers from 
instabilities [27], [28|. In atomic systems this will most 
likely lead to a phase separation into a BCS-like phase 
with equal densities and an unpaired phase with unequal 
densities. In principle, something similar could happen 
in quark matter as well [l2|, [29|, [3(| . However, because of 
long-range Coulomb forces, the different phase domains 
cannot grow arbitrarily large, and it is therefore unclear 
whether a mixed phase can exist at all. 

Instead, it seems reasonable that the matter becomes 
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inhomogeneous already on a microscopic scale by the for- 
mation of "crystalline condensates" (see Ref. [3l[ for a 
dedicated review) . The basic idea is to form Cooper pairs 
with non-zero total momentum. This has the obvious ad- 
vantage that the fermions in the pair no longer have to 
have opposite momenta, and therefore each of them can 
stay on its respective Fermi surface. In the context of 
color superconductors, this possibility has been investi- 
gated first in Ref. [HJ . The authors restricted themselves 
to a two-flavor model with a plane-wave ansatz for the 
gap function, like in the original work by Fulde and Fer- 
rell Indeed, as was shown in Ref. [33|, one of the 
instabilities which occur in gapless two-flavor color su- 
perconductors could be related to an instability against 
the formation of a Fulde-Ferrell (FF) -like condensate. 

On the other hand, since in the FF ansatz the to- 
tal momentum of the pair is restricted to a non-zero 
but constant value 2q, this pairing pattern is strongly 
disfavored by phase space in most cases. Several au- 
thors have therefore extended the ansatz to multiple 
plane waves, studying both, two- and three-flavor sys- 
tems [34], |35|, |36|, |37| • As expected, the resulting solutions 
were found to be strongly favored against the FF phase. 
However, these analyses were restricted to a Ginzburg- 
Landau approximation. This turned out to be especially 
problematic in the two-flavor case, where the Ginzburg- 
Landau functional for the energetically favored solution 
was not bounded from below [34]]. Moreover, the crystal 
structures considered so far have been restricted to super- 
positions of a finite number of plane waves whose wave 
vectors all have the same length, whereas one should also 
allow for the superposition of different wave lengths. 

The aim of this paper is to overcome the restric- 
tion to the Ginzburg-Landau approximation and to ap- 
proach the mean-field problem explicitly. The thermody- 
namic potential is then always bounded from below and 
a proper treatment leads to new insights which could not 
be obtained in the previous investigations. As a first step 
we will focus on two-flavor pairing allowing an arbitrary 
real gap functions with general one-dimensional periodic 
structure. 

For inhomogeneous ground states the mean-field prob- 
lem is already non-trivial and requires to solve the 
Bogoliubov-de Gennes equations [3g|. Only for 1 + 1- 
dimensional systems there is a good understanding of the 
mean-field ground state and the thermodynamic proper- 
ties of the system [U, HJ SI E 0, H SI- This is, 
however, lacking for higher dimensional systems and at- 
tempts have been made to simplify these equations, e.g., 
by integrating out short-range fluctuations [46|, |47| • We 
will not pursue such a direction here, but instead present 
a numerical approach to solve the Bogoliubov-de Gennes 
equations in a convenient basis. The presentation and 
derivation is elementary so that no prior knowledge of in- 
homogeneous phases is required. It turns out that at least 
for relativistic systems the regularization of the theory 
has to be addressed carefully in order to avoid undesired 
artifacts. 



The paper is organized as follows: In section |TT] we 
introduce the model we aim to investigate in a certain 
approximation and derive an expression for the thermo- 
dynamic potential in an inhomogeneous phase together 
with the corresponding gap equation. Because of its im- 
portance for inhomogeneous phases in 3 + 1-dimensions, 
we also discuss a suitable regularization scheme. In sec- 
tion IIIII we then present numerical results for inhomo- 
geneous phases with one-dimensional inhomogeneity in 
3 + 1-dimcnsions. As a prelude we discuss the homoge- 
nous (color-)superconducting and the Fulde-Ferrell phase 
first, before confronting them with results for a general 
pairing pattern. For the latter we continue by exploring 
the quasi-particle spectrum and a comparison to ana- 
lytical results obtained in 1 + 1-dimensions. Finally we 
summarize our results in section IIVI and give an outlook 
for possible further investigations. 



II. FORMALISM 

In this section we develop the general framework for 
the description of inhomogeneous color-superconducting 
phases. 



A. Model Lagrangian 

We consider an NJL-type Lagrangian for massless 
quarks q with three flavor and three color degrees of free- 
dom, 

C = q(iif) + fi)q + £ int . (1) 

We have introduced the notation tfi = /i7°, where /j is 
the chemical potential. To be precise, y, is a diagonal ma- 
trix in color-flavor space, allowing for different chemical 
potentials for different colors or flavors. 
The interaction term is given by 

C lnt =H 2^ {qi~f5T A \ A , q c ){qc ilb^A^A' q) ■ (2) 
A,A'=2,5J 

Here H is a dimensionful coupling constant and qc(x) = 
Cq T (x), where C = ij 2 "f° is the matrix of charge conju- 
gation, ta and \a> denote the antisymmetric Gell-Mann 
matrices acting in flavor space and color space, respec- 
tively. Thus, Li n t corresponds to a quark-quark interac- 
tion in the scalar flavor-antitriplet color antitriplct chan- 
nel. 

The above Lagrangian should be viewed as a typical 
example which allows for the most important pairing pat- 
terns in color superconductivity, like the two-flavor su- 
perconducting (2SC) phase and the color-flavor locked 
(CFL) phase. However, the formalism we are going to 
develop in this section is by no means restricted to this 
model. In particular we may add mass terms, and the in- 
clusion of other interaction channels is straight forward. 
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Applying standard bosonization techniques, the inter- 
action term, Eq. @, can equivalently be rewritten as 



Cint = ~ £ | {qi5TA>^A> 



qc) VAA< 



A,A' 



-<Paa< (ic 15TA Aa' q) - 2H <Paa> <PAA'} , 

(3) 

with the auxiliary complex boson fields <Paa'{ x )i which, 
by the equations of motion, 

(fiAA'(x) ^~2H q c (x)j 5 T A XA' q(x) , 

<PaA'(. x ) = 2H q(x) j 5 t a \a> qc{x) , (4) 

can be identified with scalar diquarks. 

In mean field approximation we replace these quantum 
fields by their expectation values 



( ifiAA' (x)) = A A (x) 5 A A> , 
{<Paa'( x )) = A* a (x)S A a> , 



(5) 



where the "gap function" Aa(x) is now a classical field. 
Here we assume that the condensation takes place only in 
the diagonal flavor-color components of the gap matrix, 
A = A' , as in the standard ansatz for the CFL or the 2SC 
phase. Note, however, that we retain the full space-time 
dependence of the field. 

Introducing Nambu-Gor'kov bispinors, 

(6) 

we obtain the effective mean-field Lagrangian 
C MF {x) = *(*) S~ l (x) - £ \&a(x)\ 2 , (7) 



(8) 



with the inverse dressed quark propagator 
q -l M _ ( l/i Hx)l5 

[x) i-A-(x) 75 

Here we used the more compact notation 

A(x) = £A a (x)taAa, (9) 

A 

i.e., A(x) is a matrix in color and flavor space. 

B. Thermodynamic potential 

We now consider a static crystalline structure with a 
unit cell spanned by three linearly independent vectors 
fix, a.2-, and fi$. This means, the gap matrix A (as) is time 
independent and periodic in space, 

A(x) = A(x) = A(x + fii), i= 1,2,3. (10) 



Hence, A can be decomposed into a discrete set of Fourier 
components, 



A0r)=£A 9 



-tqk-x 



(11) 



where the allowed momenta are given by the conditions 



qu 



q k -3i = 2ir N ki 



(12) 



for € 7L. These momenta form a reciprocal lattice 
(R.L.) in momentum space. 

For the bispinors ^> and \& we consider a finite quan- 
tization volume V with periodic boundary conditions. 
Working at finite temperature T and employing Matsub- 
ara formalism the (imaginary) time variable is restricted 
to a finite interval as well, < t = it < 1/T. Hence, 
the allowed energies and three-momenta both are dis- 
crete and we have the Fourier decompositions 



(13) 



Here we have explicitly taken out a normalization fac- 
tor X/s/V to have dimensionless Fourier components ^ Prl 
and ^ Pn . 

For a consistent description of the crystal, the quanti- 
zation volume should contain an integer number of unit 
cells. Without loss of generality we therefore assume that 
V is spanned by the vectors N a, where iV is a positive 
integer. Then the allowed momenta are given by 



Pn 



Pn 



p n .a l = 2n^, (14) 



with uj Pn being fermionic Matsubara frequencies and 
N n i £ TL. Comparing this with Eq. (|12|). we see that the 
three-momenta form a mesh which, in each direction, is 
N times finer than the reciprocal lattice of the crystal. 
Later we will take the infinite volume limit, N — * oo, 
where the set of allowed three-momenta becomes contin- 
uous. 

The thermodynamic potential per volume is given by 



n(r,/i) 



T 
V 



ln£(T,/z), 



where 



Z(T,/i) = V^W e 



(15) 



(16) 



is the grand canonical partition function with the Eu- 
clidean action 



S = 



,1/T , 

/ dr d 3 x £{x° = -ir, x) . 
Jo Jv 



(17) 



4 



Inserting the Fourier decompositions, Eqs. (jlip and 
(fT5)l . into Eq. and turning out the integrals, we obtain 
in mean-field approximation 



(18) 



A qk 



where 



(A. 



y a* 



75 <V, 



(19) 



is the (p m ,p n )-component of the inverse quark propa- 
gator in momentum representation. Note that in gene- 
ral S^ 1 is not diagonal in momentum space because the 
condensates A Qk couple different momenta. Physically, 
this corresponds to processes like the absorption of a 
quark with momentum p n by the condensate together 
with the emission of an antiquark or a hole with momen- 
tum p m = p n + qk ■ This is only possible because the inho- 
mogeneous diquark condensates carry momentum. In the 
homogeneous case, A(x) = const., only the momentum 
component qk = exists, and the in- and outgoing quark 
momenta are equal. While this is no longer true for our 
inhomogeneous ansatz, the fact that we consider a static 
solution still guarantees that the energy of the quark is 
conserved, see Eq. (|T2|) . This means, S' 1 is still diagonal 
in the Matsubara frequencies uj Pn . 

Since the action, Eq. (fT8"|) , is bilinear in the fields (plus 
a field independent term) the mean field thermodynamic 
potential is readily evaluated. We obtain 



with 



f> (7» 



Tr In — 

2V T 



(20) 



(21) 



where the trace is to be taken over the Nambu-Gor'kov, 
Dirac, color, flavor, and momentum components of the in- 
verse propagator. The factor ^ in front corrects for over- 
counting due to the artificial doubling of the degrees of 
freedom in Nambu-Gor'kov formalism. 

As pointed out above, the inverse propagator is diago- 
nal in the energy components. This allows us to perform 
the energy trace, i.e., the Matsubara sum in the usual 
way. To that end we write 



S pLp* = 7° {iu Pn 



T^i> m ,pJ^uj Pm .,uj Pn , (22) 



with the effective Hamilton operator 



Hp m .f>n 



{l°fn - M) 5p m ,p 



(23) 



which does not depend on io Vn . Here we have introduced 
the notation f = 7 ■ p ■ 

Since Ti is hermitian, it can in principle be diagonal- 
ized. We can then employ the formula 



T£m 



iuj Pn + Ex 
T 



f + T.n(! 



-E X /T 



(24) 



to turn out the Matsubara sum. In this way we obtain 



n (T,/i) 



-E 



(e x + 2T In (l + 2e- Ex/T 



where the sum is over all eigenvalues Ex of Ti in Nambu- 
Gor'kov, Dirac, color, flavor, and three-momentum space. 

Eq. (|25p is formally the same as for homogeneous con- 
densates. In practice, since Ti is not diagonal in three- 
momentum space, its diagonalization is of course much 
more difficult in the inhomogeneous case. However, as 
a consequence of the periodicity of the crystal, Ti can 
be brought into block diagonal form. As obvious from 
Eq. (|23p , only those quark momenta p m and p n are cou- 
pled which differ by a momentum q n belonging to the 
R.L. of the crystal. On the other hand we have seen ear- 
lier that, for a quantization volume V containing N 3 unit 
cells, the mesh of allowed quark momenta is iV 3 times 
finer than the R.L., cf. Eqs. (fl2 ]l and (fT3 )) , Therefore, Ti 
can be decomposed into N 3 independent blocks in mo- 
mentum space. The sum over the eigenstates A in Eq. (|25p 
thus separates into a sum over the different blocks times 
a sum over the eigenstates of each block. The reader 
may recognize that this structure is deeply related to 
the Bloch theorem which basically says that eigenfunc- 
tions can be labelled by a vector in the Brillouin zone 
(B.Z.) and that eigenfunctions for different vectors are 
orthogonal. 

More precisely, we write 



Pn 



Pr, 



(26) 



where q m and q n are elements the R.L. and k m and k n 
belong to the B.Z. Then p m and p n are coupled only if 
k m = k n . For each vector k n in the B.Z. we therefore 
define a projector 



(*). - 



E 

<J m ,<?„G-R.L. 



(27) 



which commutes with Ti and which projects out the block 
of coupled momenta related to k n . The effective Hamilton 
operator Ti, Eq. (f2"3")l . can thus be written as a direct sum 



(28) 



k n £B.Z. 



where 



(l°fn ~ /i) VmA ~\ m ~p n 7°75 



ffm,Pn y A; n _ pm 7075 (7Vn+M)% m ,p„ 



(29) 
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is the non-trivial part of Ti.. Here p m and p n are re- 
stricted to the corresponding subspace. 

Accordingly, we obtain for the thermodynamic poten- 
tial 



n (T,/i) = 



2T In 1 + e 



-E x (k n )/T 



k n eB.Z. 



(30) 



where E\{k n ) are the non-trivial eigenvalues of 7i(fc„ 
Finally, we can take the infinite volume limit, 



1 

V 



E 

k n eB.Z. 



B.Z. 



d 3 k 

(27T) 5 



(31) 



We then obtain 



Q (T,m) 



l/(^E{^(«) + ^>n(l + e-^)}. 

(32) 



In particular for T = 0, we have 



n (r = o lM ) = - i _/ ^E|^( fc ) 



(33) 



B.Z. 



These formulas are of course consistent with the homo- 
geneous case. In this limit the "reciprocal lattice" only 
consists of the point q = and the B.Z. is the entire 
three-momentum space. 

Before closing this section, let us give an interpreta- 
tion of the momenta in the equations above. Inspecting 
the upper right Nambu-Gor'kov component in Eq. (|2U)) . 
we see that A. Pm - Pn couples an incoming hole with mo- 
mentum p n to an outgoing particle with momentum p m . 
This means, the condensate contains a fermion pair with 
momenta — p n andp m , respectively. Writing p rn = k + q m 
and — p n = —k — q n with k € B.Z. and q m , q n G R.L., we 
see that p m — p n is just the total momentum of the pair, 
whereas 2k is the relative momentum modulo momenta 
of the R.L. 



C. Regularization 

The above expressions for the thermodynamic poten- 
tial are quartically divergent if the integral and the sum 
are left unconstrained. Therefore, we have to specify 
a regularization procedure to get a well defined result. 
Since later we want to compare the free energies of in- 
homogeneous and homogeneous solutions, it is of course 
crucial to regularize both cases in a consistent way. 

We do not want to attach a physical meaning to the 
regularization scheme. Instead we think of a local theory, 



which - at a given order in some power-counting scheme 
- can be "renormalized" by adding a finite number of lo- 
cal operators. The corresponding counter terms should be 
determinable in the homogeneous phase and be express- 
ible through physical observables. As a consequence, the 
Ginzburg-Landau coefficients, which can be derived from 
the mean-field thermodynamic potential, should depend 
on the regularization only indirectly through physical ob- 
servables, like the BCS gap at a given chemical potential. 

Given these constraints, it is not obvious how to gener- 
alize a three-momentum cutoff regularization to inhomo- 
geneous phases and as discussed in Appendix the most 
naive approach of restricting the momenta k of external 
and internal quarks to \k\ < k ot kp — h. < \k\ < kp + A 
does not meet our requirements. 

We therefore suggest a Pauli-Villars-like regularization 
scheme, which we introduce via a proper-time regulariza- 
tion of the functional logarithm in the thermodynamic 
potential (see Eq. l2"Tj) and which does therefore not rely 
on homogeneous ground states. 

In a first step we go back to Eqs. (f2Tj) and ([22]) and 
combine positive and negative Matsubara frequencies to 
get 



T^ln(z 



H) 



T 



hx(iwn - TL) {-iu> n - H) 

(34) 



where A — uj n + Ti. 2 is a hermitian operator with pos- 
itive eigenvalues. Next we replace the logarithm by its 
Schwinger proper-time representation, 

poo r l T 

In A - - / —f(r)e- TA , (35) 

JO T 

where we introduced a blocking function /(r) as a regu- 
lator. Thus, our regularization scheme is defined by spec- 
ifying f(r). 

The most simple prescription would be to put a lower 
bound in the proper-time variable, f(r) = 6(t — 1/A 2 ). 
However, as we would like to keep a structure which al- 
lows us to perform the Matsubara sum analytically, we 
prefer the function 

/(r) = l-2e- rA2 +e~ 2TA2 . (36) 

Inserting this into Eq. (|35[) . this amounts to the replace- 
ment 

In A -> In A - 2 In (A + A 2 ) + In (A + 2A 2 ) , (37) 

and we can carry out the Matsubara sum as before. Then 
the regularized version of Eq. (|32|) reads: 

1 f d 3 k 
"4 



n (T, M ) 



B.Z. 



(27T) 



E 



2 

E 

3=0 



Cj {E XJ (k) + 2T\n (l + e- E ^» T )}, 

(38) 
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where cq = Ci = 1, c± = —2, and 

E x , j (k) = y/El(k)+jA* . (39) 

Eqs. p7|) and ([38]) are reminiscent of Pauli-Villars regu- 
larization. Note, however, that according to Eq. (|39|) we 
replace the free energies in a Pauli-Villars-likc manner, 
which is not exactly the same as introducing regulator 
particles with large masses, as in the standard Pauli- 
Villars regularization. 

The two regulator terms generated by the blocking 
function Eq. (|36[) regularize quadratic divergencies. This 
is sufficient to get finite results for derivatives of the 
thermodynamic potential, like the quark number den- 
sity or the derivatives ^ n , which appear in the 

gap equations. On the other hand, since the unregular- 
ized thermodynamic potential is quartically divergent, it 
remains logarithmically divergent, even after the regu- 
larization. It should be kept in mind, however, that the 
thermodynamic potential is physically meaningful only 
up to a constant. We can therefore subtract the remain- 
ing divergency by calculating the difference to some ref- 
erence point, like the ground state in vacuum or simply 
the normal conducting phase at some given chemical po- 
tential. 



block H(k), related to the offset momentum k € B.Z., 
cf. Eq. (|29p. However, for the sake of notational brevity, 
we will drop the argument (k) in the following. 

In order to reduce the complexity, and guided by high- 
density effective theory, we want to get rid of the Dirac 
structure in Ti. We first note that the eigenvalue spec- 
trum does not change under unitary transformations. 
Choosing 

U=(l^ 5 )=-Ui = -U-\ (40) 
we obtain W = WHU with 

v i = ( (l°fn ~ M) 8? m ,p n A Pm -p„ \ 

(41) 

In homogeneous phases, a standard method to diagonal- 
ize the remaining Dirac structure is to employ energy 
projectors, 

A± = i(l± 7 °?5), (42) 
where p = p/p with p — \p\ . We can then rccxprcss 

ffn ~H= (Pn - M)At - (Pn + , ( 4 3) 



D. Simplified model 

The general framework for the mean-field thermody- 
namic potential derived in Sec. IIIB| together with the 
regularization procedure suggested in Sec. Ill C| is one of 
our central results. It may serve as a starting point for 
extensive studies of the phase diagram of strongly in- 
teracting matter in future investigations. In most of the 
remainder of the present article, we will illustrate the 
power of our approach by numerical examples in a sim- 
plified version of our model. 

1. Dirac structure 

In a first step, we derive an approximation to our 
model, which should be valid at high densities. Start- 
ing point is the effective Hamilton operator H, Eq. (|23p. 
As we have discussed, Ti is block diagonal in momen- 
tum space, and we can in fact concentrate on a single 



and Ti' can be decomposed into a positive and a negative 
energy part, which act on orthogonal subspaces of the 
Hilbert space. One can therefore find a new basis, where 
these positive and negative energy parts decouple. 

However, as the are only projectors for states with 
the same momentum direction p, the method described 
above, does not work exactly in inhomogeneous phases, 
where different momenta are coupled by the gap func- 
tions. In fact, we can still remove the Dirac structure 
from the diagonal momentum components in this way, 
but at the same time new Dirac structures appear in the 
off-diagonal components. 

We can work this out explicitly by performing a uni- 
tary transformation Ti" = V'TC'V which diagonalizcs 
the Dirac part of the diagonal momentum components 
Tt'fi m s m ■ The result reads 

H" = H" + Hz , (44) 

where 



{Ti 1 )p m ,p n 



A* 

V 



-(Pr, 



■(jPm + A 1 ) ( V„ 

A! 



(45) 



Pm —Pn 

(Pm + M) $p m ,p n / 
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is diagonal in Dirac space. Here we have reordered the 
lines and columns to make the diagonal structure ob- 
vious. The non-diagonal 2x2 blocks now describe the 
Nambu-Gor'kov structure. The 2x2 identity matrix on 



the right is related to the spin degeneracy of the problem 
and is, thus, part of the Dirac structure. 

In this basis, the remaining part of the transformed 
Hamiltonian is given by 



(^2)p m ,Pn 



/ 
A 



Pn -Pr, 







—pn 



-A? 



Pm -~Pn 



— p n 



V- A : 



Pn-P„ 



A 



Pn-Pr, 



■Pn — 1) 12 + i{ftm X Pn) ■ <?) , (46) 



/ 



where cr denotes the Pauli matrices, indicating a non- 
trivial spin structure. Obviously, TL' 2 ' is not diagonal in 
Dirac space. However, the components of H'2 vanish for 
parallel momenta, p m = p n . This reflects the fact that 
states with parallel momenta have the same energy pro- 
jectors . 

In the following, we will neglect Ti'^- We expect that 
this is a good approximation at very high densities where 
the physics is dominated by collinear scattering near the 
Fermi surface, \p m — p n \ -C p m ~ p n ~ p>- Neglecting the 
antiparticle contributions as well, we obtain 

Sio(I» 

=i/|E(«^»('«-*)) 

B.Z. A 

+ regulator terms , (47) 

where E\(k) are now the eigenvalues of the "high-density 
effective Hamiltonian" 

(48) 

for a given offset momentum k. 



2. Color-flavor structure 

The second simplification concerns the color-flavor 
structure of the model. In Sec. Ill Al we have chosen a form 
of the gap matrix which is capable to describe the most 
common phases in color superconductivity, see Eq. 
For instance, the CFL phase corresponds to the case 
A 2 = A 5 = A 7 = A. 

In the following, we restrict ourselves to the 2SC pair- 
ing pattern, which is defined by A2 = A and A5 = A7 = 
0. Moreover, since we are only interested in the effect of 
the pairing relative to the normal conducting phase, we 
can omit those colors and flavors which do not participate 
in the pairing (i.e., strange quarks and, using standard 



nomenclature, blue quarks). Thus, the remaining Hamil- 
tonian is a 4 x 4 matrix in color-flavor space which can 
trivially be decomposed into four separate blocks. 

We assume that the chemical potential may be diffe- 
rent for up and down quarks, 

Hu = p + Sfi , fi d = p- 5/j, , (49) 

but does not depend on color. We then obtain 

'Hhde = / Ha,Su © Ti—AJn ffi Ha,-j^ © "H-a-Su , (50) 

where 

("HA,Sfj,)p m ,p n = 

( (Pm -fi-Sfji) Sf m> ff n &p m -p n \ 

V A * Pn -p m " {Pm ~ P + Sfi) <5p m ,p„ J 

(51) 

Of course, the eigenvalues of Ha,5u do not depend on 
the overall sign of A. Moreover, at least in all cases to be 
considered in this article, replacing <5/i by — <5/i amounts 
to a replacement of the eigenvalue spectrum {E\(k)} by 
{— E\(k)} (see Appendix |B|) . Hence, since the thermo- 
dynamic potential depends only on the moduli of the 
eigenvalues, each of the four blocks contributes equally, 
and we only need to determine the eigenvalues of TLa,s^,- 



3. Gap equations 

Summarizing our main equations, including the ap- 
proximations introduced above, the regularized mean- 
field thermodynamic potential is given by 

n(T, p, Sti) = n (T, p, sn) + ± I 2 > ( 52 ) 
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with 



noCT,fcffA0 = -2 /(^aE 



B.Z. 



Y J C ] {E Kl (k) + 2T\n (l 



-E X:] (k)/T\ 



3=0 



where cq = Ci = 1, 
+ J A 2 as in E q. 

values of Wa,5/j(^), Eq 
main focus, this becomes 



(53) 



n (T = 0,/2,^) 



ci = —2, and E\ tJ (k) — 
15]). E\{k) are now the eigen- 
H). At T = 0, which will be our 



B.Z. 



(2 



A j=0 



(54) 



As discussed earlier, in order to get a finite result for 
f2, one still has to subtract an infinite constant. In this 
paper, we will always consider the free-energy difference 
to the normal phase, i.e., the phase with A = at the 
same value of T, ft, and dfi. This quantity is finite. 

Finally, before coming to the numerical results, we 
want to discuss the gap equations. As in the homogeneous 
case we have to minimize the thermodynamic potential, 
which means that we have to solve the equations 







dfl 



dn 



^Qk 

AH 



(55) 



for all Fourier components. From Eq. (|53[) . we obtain 

dn 



dA 



Qk 



B.Z. 



d 3 k t-, dE x 
(2tt) 3 ^ <9A* 

v 1 A 1 k 



^c j (l-2n(Ex, j ) 

3=0 



E 



Ex 



, (56) 



where n(E) = [exp(E/T) + is a Fermi function. To 
evaluate this further, we use that Ex are the eigenvalues 
of Ha,5^- This means, there are unitary matrices U so 
that 

Ha,5„ = U- 1 Ha,s i .U (57) 
is a diagonal matrix with (Ha.s^xx = E\. Hence, 
dEx d 



dA 



Qk 



dA 



Qk 



where P qk = 



OA' Wa,5a» 

Ik 



(58) 

is a known A- and k- 



independent matrix (see Eq. (|6ip below). Note that the 
terms related to the derivatives of and U cancel each 



other. Writing the matrix U = (wi, . . . ) in terms of the 
eigenvectors wx to the eigenvalues Ex, this yields 



dEx 
dA* 



(59) 



In fact, this formula is nothing but first-order perturba- 
tion theory for the modification of the eigenvalues under 
a small perturbation STLa.&h = Pq^^-Qk - If we insert this 
into Eq. (l56|) . we can then rewrite Eq. ([55]) into the gap 
equations 



B.Z. 



q k W X 



J2c j (l-2n(E x , j j)^L. (60) 



These equations are the basis for our numerical analysis. 

At this point, we would like to note that the matrix P qk 
connects momenta differing by qk- In fact, from Eq. (|5ip 
we get 



{P qk )i 





1 



(61) 



Thus, denoting the Nambu-Gor'kov components of the 
eigenvectors as wx — (ux, vx) and indicating the momen- 
tum components explicitly, the gap equations read 

/d 3 k / \ 

7^3 E K)pW* («Ak (l - MEx)) , 

(62) 

where we have omitted the regulator terms for simplicity. 
Fourier transforming this convolution with the conven- 
tions given in Eqs. (I11I13P we get 



B.Z. X 



x) ux{x)(l-2n{Ex) 



where V is the volume of the unit cell. This is pre- 
cisely the selfconsistency condition for the Bogoliubov- 
de Gennes equation [38| when again exploiting Bloch's 
theorem. It is obvious that this relation is not restricted 
to the simplified model we discussed here, but an exten- 
sion to the general case is straight forward. Note however 
that our regularization procedure and therefore also our 
numerical calculations are tied in momentum space. 



III. NUMERICAL RESULTS 

In this section we want to discuss numerical calcula- 
tions performed within the simplified model of Sec. Ill Dl 
We restrict ourselves to T — and to a fixed average 
chemical potential p, = 400 MeV. Then S/i is the only re- 
maining external variable and we will drop the arguments 
T and \x in the following. 
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Our model has two parameters, namely the coupling 
constant H and the cutoff parameter A. We remind that 
A restricts the free energies and not the momenta. Thus, 
the most relevant excitations around the Fermi surface 
are always included, and there is no need for A to be 
larger than the chemical potential. Of course, A should 
be considerably larger than the gap. Having fixed the 
cutoff, we will express the coupling constant H through 
the corresponding value of the BCS gap. 

For given model parameters and Sfi, the thermody- 
namic potential, as defined above, depends on the gap 
function A. Our main goal is to find the most stable so- 
lution, i.e., the minimum of f2 with respect to A. At a 
given periodicity of the crystal, this corresponds to min- 
imizing Q with respect to the Fourier components A qk , 
i.e., to finding the most favored solution of the coupled 
set of gap equations, Eq. (|60|) . In addition, we should 
vary the periodic structure itself, i.e., the basis vectors of 
the reciprocal lattice. Obviously, this is a very involved 
problem, which is beyond the scope of the present paper. 

Therefore, as a first step, we restrict ourselves to one- 
dimensional crystalline structures, i.e., to gap functions 
which vary periodically in one spatial direction q, but 
stay constant in the two spatial directions perpendicular 
to q. Moreover, we only consider real gap functions. In 
spite of these restrictions, we find an interesting class 
of solutions, which will be discussed in Sec. IIII CI and 
thereafter. First, however, we briefly discuss the BCS and 
the Fulde-Ferrell solutions within our model in order to 
provide a basis and to make contact to the literature. 

A. BCS phase 

The BCS solution corresponds to the limit of spatially 
homogeneous condensates. In this case the Brillouin zone 
is the entire three-dimensional space and we have 

A gfc = AS qk<0 . (64) 

Therefore, the effective Hamiltonian Eq. (f5~lj) is diagonal 
in momentum space, and we obtain the eigenvalues 

\E ± {k)\ = W(k-fi) 2 + \A\ 2 ±S^\ . (65) 

Inserting this into Eqs. ([52|) and (|54|) we recover the well- 
known result that for Sfi < |A| the unregularized part of 
the thermodynamic potential, 

r f pi. |A|2 

Slunre 9 (Sfl) = -4 J — ^ y/ (k W + I A| 2 + LL , 

(66) 

is independent of Sfi When the Pauli-Villars reg- 
ulators are included, this does no longer hold exactly. 
As a consequence, the value of A which solves the gap 
equation dft/dA* = is weakly Sfi dependent. For this 
reason we define Abcs to be the gap at Sfi = 0. For 
Abcs = 80 MeV and A = 400 MeV, which will be our 
standard choice of parameters, we then find that A in- 
creases by about five percent, when Sfi is varied between 



and 0.8 Abcs- (Higher values of Sfi are irrelevant for 
the BCS phase, as we will see below.) For smaller gaps 
or larger values of the cutoff the effect is even smaller. 

When Sfi is increased, the BCS phase eventually be- 
comes unfavored against the normal phase. The corre- 
sponding phase transition is first order. In the weak- 
coupling limit it occurs at Sfi = Abcs/v^, as shown 
already in 1962 by Chandrasekhar and Clogston [22|. For 
stronger couplings, the BCS phase can sustain somewhat 
larger values of Sfi. We will come back to this in the fol- 
lowing subsection. 

B. Fulde-Ferrell solutions 

The FF phase corresponds to the case that the gap 
function is a single plane wave, 

A{x) = Ae 2tff . (67) 

This means, the B.Z. is infinite in the directions per- 
pendicular to q, but finite in g-direction with length 2\q\. 
The momentum space representation of the gap matrix 
is given by 

A 9fc =A%, 29 -. (68) 

Inserting this into Eq. (|5T|) one finds that the effective 
Hamiltonian is still block diagonal, however with shifted 
blocks where two momenta are coupled. These blocks can 
be written in the form 

( t+ -r J "- t _ + v {( ,). c») 

with k± = \k±q\. Moreover, the sum over all blocks sim- 
ply amounts to extending the /c-integration to the entire 
space. Hence, the thermodynamic potential reads 

"CM) = -*J 7^3 £ (1^(^)1 + reg.) + ^ , 

(70) 

where 

E ± {k) = k -±^- ±\/(^^-m) 2 + |A| 2 -^ 

(71) 

are the eigenvalues of the matrix (f69|) . Apart from the 
Pauli-Villars regulators in Eq. ([TO")) , this is again a stan- 
dard result. 

In order to find the most favored FF solution, we must 
minimize the thermodynamic potential with respect to 
|A| and Typically, one finds that |<f| is of the or- 
der of 0.9 Abcs while |A| is considerably smaller than 
Abcs- 1 Numerical examples for Abcs — 80 MeV and 



In the weak-coupling limit, \q\ = 0.906 Ages and |A| = 
0.23 Ages at the Chandrasekhar-Clogston point [lj,0|- 
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A = 400 MeV will be discussed in the context of Figs. [5] 
and [H 

In competition with the normal phase and the BCS 
phase, this optimal FF solution is favored only in a small 
window in Sfi. This is shown in Fig. [T] where the phase 
boundaries are displayed as functions of Abcs for two 
different values of the cutoff. The phase transitions are 
of first order for the BCS-FF transition (dashed lines) 
and of second order for the FF-normal transition (solid 
lines). Indeed, in the weak-coupling limit, one expects a 
first-order phase transition from the BCS phase to the 
FF phase near the Chandrasekhar-Clogston value <5/i = 
Abcs/V^ followed by the second-order phase transition 
to the normal phase at 5fx = 0.754 Abcs 0, 0j- These 
values are indicated in the figure by the thin horizontal 
lines. Obviously, the model results come close to these 
limits for Abcs ~ > 0, whereas for stronger interactions 
we again find deviations. 

The BCS-FF phase boundaries are almost identical 
with the corresponding BCS-normal boundaries because 
the free energy difference between BCS phase and normal 
phase is a much steeper function of Sfi than the free en- 
ergy difference between FF phase and normal phase (see 
Fig. [4J. This is a standard result, which was also found, 
e.g., in Ref. [13] where a different regularization scheme 
was used. 

The FF-normal phase boundary, on the other hand, 
behaves rather differently from the result of Ref. [32| . 
There, it was found that the critical 8fi decreases with in- 
creasing Abcs- As a consequence, the BCS-FF boundary 
and the FF-normal boundary intersect at some Abcs ~ 
100 MeV, and there is no stable FF solution at stronger 
couplings. In contrast to this, we find that the FF-normal 
phase boundary runs almost parallel to the BCS-FF 
phase boundary so that the width of the FF window stays 
approximately constant. This means, as long as we do not 
consider other inhomogeneous phases, the regime of sta- 
ble FF solutions extends to very large values of A bcs m 
our regularization scheme. 



General solutions for real one-dimensional gap 
functions 



We now turn to the main part of this work, being the 
analysis of ground states when allowing for a real one- 
dimensional gap function as the order parameter. For 
simplicity, we will often call these solutions "general so- 
lutions" , although there should be, of course, even more 
general solutions with higher-dimensional or complex gap 
functions. 

In the first step we limit ourself to periodic gap func- 
tions of the form 



0.82 
0.8 
0.78 



FF-nor, A= 400MeV 
FF-nor, A= 8()0MeV 
BCS-FF, A= 400MeV 
BCS-FF, A= 800MeV 
8n= 0.707 A BCS 
8|X= 0.754 A 




60 80 
[MeV] 



120 



FIG. 1: Stability window in Sfi for the FF phase in compe- 
tition with the BCS phase and the normal phase as a func- 
tion of the coupling, parameterized by Abcs- The BCS phase 
is favored below the dashed lines, the normal phase is fa- 
vored above the solid lines. Results for two different cutoffs are 
shown. The thin horizontal lines indicate the weak-coupling 
limits. 



components A^^ and afterwards with respect to \q\. 

Comparing Eq. (|72]) with Eqs. fTl]) and (fl"2j) . we find 
that qk — 2kq. The (admittedly somewhat unnatural) 
factor of 2 was introduced to give q the same meaning 
as in the FF ansatz, Eq. ([67|) . where this factor is the 
standard convention. 

Again the B.Z. is infinite in the directions perpendic- 
ular to q and finite in the (f-direction with length 2\q\. 
Without loss of generality, we can assume that q points 
in z-direction, q = qe z . Then one period corresponds to 
z = 

q 

The restriction to gap functions which are real in coor- 
dinate space means that the Fourier components satisfy 
the relation A^ = At_ k . Moreover, without loss of gen- 
erality, we can always choose the origin to be located at 
a maximum of the gap function. A(z) is then an even 2 
function, and the Fourier components are real. 

In the following we fix the model parameters to 
A = 400 MeV and a coupling strength corresponding 
to Abcs — 80 MeV. As before, we consider T = and 
fl = 400 MeV. 

In Fig. [2] we present examples of one-dimensional gap 
functions A(z) we obtained by minimizing the thermody- 
namic potential at <5/i = 0.7 Abcs for different fixed peri- 
ods. For convenience, z is measured in units of -nj Abcs 
so that one period is given by (q/ Abcs)^ 1 ■ At q ~ Abcs 



A(x) 



k 



Aff k e 



likq-x 



(72) 



For each period given through |g| we will then minimize 
the thermodynamic potential with respect to the Fourier 



2 In our numerical analysis, we only found solutions which are 
symmetric under reflections at a plane perpendicular to q going 
through a maximum or minimum. So far, we cannot exclude that 
other solutions exist which do not have this symmetry. In this 
case A(z) would of course not be an even function. 
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FIG. 2: The gap function in coordinate space at Sfi = 
0.7Abcs for different fixed values of q. 
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FIG. 3: Difference between the thermodynamic potentials of 
the general inhomogeneous phase and the normal phase as a 
function of q for different values of Sfi. 



the gap function appears to be sinusoidal. For larger pe- 
riods, however, a new feature becomes apparent: the for- 
mation of a soliton lattice. Especially for q = 0.1 Abcs, 
we see that the gap function stays nearly constant at 
zLAbcs for about one half-period and then changes its 
sign in a relatively small interval. The q — 0.2 Abcs so- 
lution behaves qualitatively similar, but has a shorter 
plateau. Remarkably, the shape of the two functions is 
almost identical in the transition region where the gap 
functions change sign. This remains even true for the 
q = 0.5 Abcs solution, which is kind of an extreme case 
with no plateau and only transition regions. We may thus 
interprete these transition regions as very weakly inter- 
acting solitons, which are almost unaffected by the pres- 
ence of the neighboring (anti-) solitons as long as they 
do not overlap. 

These features will be discussed in more detail in 
Sec. IIII El The main result is that the gap functions are 
characterized by two independent scales. The first scale, 
q, determines the period of the lattice and thereby the 
distance between the solitons. The second scale, Abcs, 
is not only the amplitude but also determines the shape 
of the solitons, which is practically independent of q. In 
fact, even for the sinusoidal solution at q = Abcs, the 
slope at the zero crossing is almost the same as for the 
q = 0.1 Abcs solution, although the period differs by one 
order of magnitude. 

For each Sfi, with the solutions for the different cho- 
sen values of q at hand, we now have to minimize the 
thermodynamic potential in q in order to determine the 
energetically preferred ground state. This is illustrated 
in Fig. [21 where the difference 5£l between the thermo- 
dynamic potential of the inhomogeneous phase and the 
normal phase is displayed as a function of q for three dif- 
ferent values of <5/i. We observe that the preferred value of 
q rises from small values towards Abcs when increasing 
(5/i. Also the shape of the potential changes. 



We want to point out that the BCS phase is included 
in this plot as q = 0. It is interesting to note that the 
inhomogeneous phase is the preferred phase already at 
5fx = 0.70 Abcs, i- e - 5 significantly below the transition 
between the BCS phase and the FF phase (cf. Fig. [I]). 
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FIG. 4: Difference between the thermodynamic potentials of 
different phases and the normal phase as functions of Sfi: BCS 
phase (solid line), general inhomogeneous phase (dashed line), 
and FF phase (dotted line). 

This is also seen in Fig. QJ where 8tt is displayed as a 
function of Sfi for the general solution in comparison with 
the BCS phase and the FF phase. The shown results for 
the inhomogeneous phases correspond to the preferred 
values of q at each 6fi. We find that the LOFF window, 
i.e., the interval of 5/i for which a general inhomogeneous 
ground state is energetically favored, has almost doubled 
compared to the FF window. This is due to the fact that 
the interval is expanded towards smaller values of <5/x, the 



12 



upper end is almost unchanged. Related to this, the free 
energy gain of the general solution is much larger than 
that of the FF phase. 

However, the most striking features which are visible 
in Fig. [4] are the orders of the phase transitions. Whereas 
for the FF solutions, the phase transition to the BCS 
phase is first order and to the normal phase is second 
order, it is just the other way around for the general 
inhomogeneous solutions: Here we find a first-order phase 
transition to the normal phase, whereas the transition to 
the BCS phase at the lower end appears to be continuous. 



100 




BCS 



FIG. 5: The energetically preferred value of q in the general 
inhomogeneous superconducting phase (solid) and for the FF 
solutions (dotted) as functions of Sfi. The solution for the FF 
phase stops at the transition to the normal phase. 



The latter is possible because our space of possible 
solutions allows for a natural connection between these 
phases via the formation of a soliton lattice. This becomes 
more clear in Fig. [5] where the energetically favored values 
of q are displayed as functions of Sfi. The solid line and 
the dotted line correspond to the general inhomogeneous 
phase and to the FF phase, respectively. In the FF phase, 
q is always of the order of Abcs- Similar values of q are 
also found for the general solutions in the upper part of 
the LOFF window, S/i > 0.75 Abcs- With decreasing Sfi, 
however, the preferred q decreases to arbitrarily small 
values and eventually goes to zero at Sfi « 0.695Abcs- 

We thus arrive at the following picture: With lower- 
ing Sfi, the period of the gap function increases and we 
eventually obtain a soliton lattice with increasing dis- 
tance between the solitons, i.e., with constant plateaus 
of increasing length, see Fig. [2] At the critical point, the 
length of the plateaus diverges, and the inhomogeneous 
phase is continuously connected to the BCS phase. Also 
notice that, although the transition is continuous, the 
slope of the function q(5fi) changes dramatically when 
q/Ascs comes to the order of 0.5. As we have seen in 
Fig. [21 this is just the regime, where the more sinusoidal 
solutions go over into a soliton lattice. 
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FIG. 6: The amplitude of the gap function A(z) at the en- 
ergetically preferred value of q at given Sfi for the general 
inhomogeneous phase (solid line) and the FF phase (dotted 
line) . 



In Fig. [5] the amplitude of the general inhomogeneous 
gap function is displayed as a function of Sfi (solid line). 
Since the transition to the BCS phase is continuous, the 
amplitude becomes equal to the BCS gap at the lower end 
of the window. 3 With increasing Sfi, the amplitude de- 
creases, but it remains of the order of Abcs m the entire 
interval. For comparison, we also show the amplitude of 
the FF solution (dotted line). In contrast to the general 
solution it is always considerably smaller than Abcs and 
goes to zero at the transition point to the normal phase. 
Thus, in agreement with our conclusions from Fig. |H the 
transition to the normal phase is of second order for the 
FF phase, but first order for the general inhomogeneous 
phase. The fact that the amplitude of the general solu- 
tion is considerably bigger than the amplitude in the FF 
phase also reflects our earlier observation that the general 
solution is energetically much more favored (see Fig. 0}. 

At first sight, a first-order phase transition from the in- 
homogeneous phase to the normal phase seems to be in 
contradiction with Ginzburg-Landau investigations. As 
already discussed in the context of Figs. [5] and [5] and as 
will be detailed in section lHI El the shape of the gap func- 
tion becomes more and more sinusoidal with increasing 
Sfi, i.e., increasing q. A sinusoidal (or so-called antipo- 
dal) gap function has been investigated in a Ginzburg- 
Landau analysis, showing that the transition from the 
inhomogeneous phase to the normal phase is of second 
order [3l|,[34]]. Since in the vicinity of a second-order tran- 
sition the Ginzburg-Landau approximation is expected to 
converge to the exact mean-field result, an explanation 



3 We remind that, because of the regularization terms, the BCS 
gap at Sfi > is slightly larger than Ages- 
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FIG. 7: Difference between the thermodynamic potential of an 
inhomogeneous phase with a sinusoidal gap function and the 
normal phase as a function of the magnitude of the gap func- 
tion. Here we have chosen q = 0.9Abcs and Sfj = 0.775Abcs ■ 



of this discrepancy is needed. 

To clarify this issue, we determine the thermodynamic 
potential for a gap function 

A(z) = A cos cos(2gz) (73) 

with fixed q, as being used in the Ginzburg-Landau anal- 
ysis. An example is shown in Fig. which corresponds 
to q = 0.9 Abcs and Sfi — 0.775 Abcs- We see that the 
thermodynamic potential has two minima in this sub- 
set of possible solutions. The more shallow one at lower 
magnitudes of A cos may be described by the Ginzburg- 
Landau analysis. However, in order to get the second min- 
imum, one needs at least to include terms of the order 
A^ os , whereas in Refs. f3~Lll3~ij at most terms of the order 
A® os were included. Moreover, investigations in one spa- 
tial dimension have revealed that a gradient expansion 
of the thermodynamic potential may even not exist at 
all (48J. Therefore it is unclear whether the second min- 
imum is determinable in a Ginzburg-Landau approach, 
even if higher powers in the order parameter are included. 

D. Quasiparticle spectrum 

Next we want to discuss the excitation spectrum of the 
quasiparticles and its impact on the shape of the energet- 
ically preferred gap functions. For a given vector k in the 
B.Z. we have a discrete eigenvalue spectrum of the as- 
sociated Hamiltonian Ha,s^W, Eq. (j51~j) . which depends 
on k via its allowed momenta, see Eqs. (I26H29[) . As the 
inhomogeneous gap functions break the rotational sym- 
metry of the system, this eigenvalue spectrum does not 
only depend on the modulus of k, but also on its direc- 
tion. Therefore, we will present the spectrum by showing 
selected slices through the B.Z. 



We note that the eigenvalue spectrum of Ha,8^ de- 
pends on Sfj in two ways: First, there is an explicit Sfj 
term in the diagonal Nambu-Gor'kov components and, 
second, there is an implicit Sfj dependence through the 
Sfj dependence of the gap functions. For fixed gap func- 
tions, the explicit Sfj, terms simply shift the eigenvalues 
of Ha,o by —Sfj,. We will therefore take this trivial effect 
out and show the eigenvalue spectrum of TIa,o, which we 

will denote by E^ (k) . This spectrum then still depends 
on S/i through the gap functions. The gap functions, in 
turn, depend on Sfj mainly through the variation of \q\ 
(see Figs. [5] and [5]). Therefore, we present the spectra at 
fixed values of \q\. As we will see below, the remaining 
Sfj dependence is very weak, i.e., the spectra are quite 
generic. 

In Fig. [8] we show two examples of the excitation spec- 
trum at a fixed value of the momentum k± perpendic- 
ular to q as a function of the momentum k z along the 
direction of q. For a real gap function, the eigenvalue 
spectrum possesses two symmetries: First, as shown in 
Appendix [B] the eigenvalues Ejj^ at given momentum 

k in the B.Z. appear in pairs (Ej?\ —E^). Second, the 
eigenvalue spectrum at k z and 2\q\ — k z coincide. This is 
related to the fact that the gap functions are even func- 
tions in z, i.e., parity is unbroken by the condensate. 4 

In the left panel of Fig. [5] we show the spectrum for 
k,± = 0. We see that there are four low- lying excitations 
with free energies well below Abcs- These modes are re- 
lated to each other by the symmetries discussed above, 
i.e., there is in fact only one non-trivial solution at low 
energies. This feature is also known from analytical in- 
vestigations in one spatial dimension and is related to 
the soliton [4^]. In addition to these modes, there are 
(infinitely many) other excitations with higher free ener- 
gies. These higher-lying modes are clearly separated from 
the low-lying ones, i.e., there is a gap in the excitation 
spectrum, at least as long as we keep k± fixed. 5 In con- 
trast to the BCS phase, the gap does, however, not start 
at vanishing free-energies. 

The right panel of Fig. [5] corresponds to k± — p, — 
400 MeV. Most features of this spectrum are the same 
as in the previous example. The main difference is the 
fact that now the low-lying excitation does not change 
its sign when k z is varied. As a consequence, the positive 
and the negative solutions do not cross and there is an 
additional gap around vanishing free energies. 



4 This statement holds for our particular choice of the coordinate 
system. More generally, the gap functions are symmetric under 
reflection at a plane perpendicular to q going through a maximum 
or minimum. See also footnote [2] 

We remind that here we are discussing the eigenvalue spectrum 
EY*' (k) of Wa,i5h=o- The true spectrum of Hhde ls obtained by 
shifting the results by 8(1 upwards and downwards, see Eq. 1501 . 
The gap may then disappear in some cases. This is, however, 
irrelevant for the discussion later on in this section, which is 
based on the existence of a gap in the spectrum of Ha o- 



14 




FIG. 8: Excitation spectrum in the B.Z. as a function of the momentum k z parallel to q for fixed perpendicular momentum 
k± = (left panel) and k± = 400 MeV (right panel). Shown are the eigenvalues E^\k) of 7Ya,«h=o for the gap functions of the 
energetically preferred ground state at fixed |q| = 0.9 Abcs and Sji = 0.7 Abcs- The spectrum E\(k) of TChde is obtained by 
shifting the lines by <5/i and — <5/i. 



The gap structure as a function of k±^ is visualized 
in Fig. [5] for two examples. The shaded areas indicate 
the free energies which can be reached by at least one 
mode at the given value of k± , when k z is varied over all 
possible values. Accordingly, the white areas correspond 
to the gapped regions. We see again that such gaps exist 
in the excitation spectra when k± is kept fixed. There is, 
however, no gap when all values of k± are considered. 

The left panel of Fig. [9] corresponds to |g| = 0.9 Abcs, 
i.e., to the same period as the examples shown in Fig. [5] 
Again, we see that the low-lying modes are restricted to 
a single band around zero free energies for lower values 
of k± , wheres at higher values of fcj_ this band splits into 
two. Comparing this with the right panel of Fig. [51 corre- 
sponding to \q\ — 0.2Abcs, we see that these features re- 
main qualitatively unchanged. However, we observe that 
the bands of the low-lying modes get squeezed consider- 
ably when going from larger to smaller values of \q\, i.e., 
when separating the solitons by stretching the lattice. In 
fact, we expect that in the limit \q\ — > and at low values 
of k± , the modes associated with the soliton are restricted 
to exactly zero free energy, whereas the continuum starts 
at \E\ = A B cs- 

A gap in the complete excitation spectrum has im- 
portant consequences for the dependence of the thermo- 
dynamic potential and the gap functions on 5/i: Con- 
sider the eigenvalues Ef\k) of the unshifted Hamilto- 
nian Ti.A,Sfj,=o(k). Then, as pointed out above, the eigen- 
values of HA,5ft are simply given by E x = E[ 0) - Sfx. 
Thus, since the E^ always come in pairs (e{°\ —E^), 
we get from Eq. (|54l 

nop/*) = -2 j |^ £ - m + + ^|) , 

B.Z. E { ° ] >0 

(74) 



where we have dropped the regularization terms for sim- 
plicity. This can be written as 

fio(ty) =^ 1) +n? ) (<5/x) , (75) 
with a Sfi independent part 

^=- 4 /(03 E 1^)1 (76) 

B.Z. E < -° ) >5n 

and a <5/i dependent part 

^ 2) ^) = - 4 /(0 E * 

B.Z. 0<Bf ) <5n 

/d 3 k 
j^N s (k) ee -4S M n s . (77) 

B.Z. 

Here N s (k) is the number of positive eigenvalues E x °\k) 
smaller than dfi. 

Now suppose, there was a complete gap in the spec- 
trum in some energy interval, so that N s (k) is constant 
for all values of k, when Sfi is varied within this gapped 
interval. If n s — 0, i.e., if Sfj, is smaller than the smallest 
positive eigenvalue, the thermodynamic potential and, 
hence, the gap function are obviously dfi independent. 
For n s > 0, on the other hand, the thermodynamic po- 
tential is <5/i dependent. However, even in this case, the 
gap functions are still Sfi independent, when dfi is varied 
within the gapped interval. This can be inferred from the 
gap equations, — 0. A S/i dependence of the gap 

function must then be due to the variation , . , because 
all other terms are Sfi independent. On the other hand, 
an infinitesimal variation of the gap function can only 
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FIG. 9: Superposition of the eigenvalue spectra (k) of TCa,S/j,=o in the B.Z. along the momentum k z parallel to q as 
functions of the perpendicular momentum k± (shaded areas). The spectra have been obtained for the gap functions of the 
energetically preferred ground state at S/j, — 0.7Abcs for fixed periods \q\ — 0.9A_bcs (left) and \q\ = 0.2A_bcs (right). 



lead to infinitesimal changes in the eigenvalue spectra. 
Therefore, if S/i lies in a gapped region, this variation 
will not change the numbers N s (k), which are, by defini- 



tion, integer numbers. Consequently, 
gap function is Sfi independent. 6 
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FIG. 10: Gap functions at S/j, — 0.69Abcs and 5/j, 
0.79Abcs , and their difference for |<j| = 0.8A_bcs- 



In our case the excitation spectrum is not completely 
gapped. Nevertheless, when varying Sfj, in the vicinity of 
0.7 Abcs, the only eigenvalues interfering are those asso- 
ciated to the solitons and with momenta |fc| > \x. There 



6 In our model, this property is slightly spoiled by the regulariza- 
tion, as already discussed for the BCS gap. 



number is however relatively small and they are outside 
the Fermi ball. We can therefore expect their influence 
on the value of the thermodynamic potential and, hence, 
on the shape of the gap function to be very small. 

To illustrate this property we present the gap func- 
tion at Su = 0.69 Ascs an d Su = 0.79 Abcs for |<f| = 
0.8 Abcs m Fig. |TD] The gap functions, though taken at 
the lower and upper end of the LOFF window, appear 
almost identical. We conclude that the shape of the gap 
functions in the LOFF window for given \q\ is almost 
independent of Sfi. 

We also note that the spatially averaged density dif- 
ference between up and down quarks is given by 



(Sn) 



( n d) 



We thus get from the above equations 



{Sn) 



5jJL 



di 



(78) 



(79) 



where we have used that the implicit 5fi dependence 
through the gap functions drop out because of the gap 
equations. Moreover, as we have seen before, n s is al- 
most independent of Sfj,, if Sfi is varied within the nearly 
gapped region. Therefore the second term on the right 
hand side is small and we find that the density difference 
is approximately proportional to the density of solitonic 
states. 



(Si 



4n, 



(80) 



In particular, this implies that, unlike in the BCS phase, 
the density difference does not vanish if low-lying soli- 
tons are present. On the other hand, (Sn) is almost Sfi 
independent when Sfj, is varied in the gapped regime. 
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This gives rise to the very intuitive picture that the 
excess quarks are sitting in the solitons, whereas in the 
BCS-like plateaus of the gap functions the densities of 
up and down quarks are nearly equal. Here we should 
keep in mind that the period of the gap function and, 
hence, the soliton density, was kept constant in the above 
discussion. This explains why we found (Sn) to be <5/x 
independent. If we do not fix the period, but minimize 
the thermodynamic potential with respect to \g\, as in 
Fig. [5l ri is no longer linear in 5/j,, and (Sn) becomes S/i 
dependent. In fact, we can turn this argument around 
and explain the Sfj, dependence of |g| by the preference 
of the system to accommodate more excess quarks with 
increasing <5/i. 

E. Comparison with analytical results in one 
spatial dimension 

It is quite instructive to compare our results for the 
gap functions in three spatial dimensions with the an- 
alytically known solutions of the mean-field problem in 
one spatial dimension. 

In 1 + 1-dimcnsions the mean-field problem can be ap- 
proached in different ways and its first solution goes back 
to Peierls [39]. In terms of inverse scattering theory, the 
selfconsistently determined gap function has to generate 
a reflectionless (for a single bound state) or finite-gap 
potential in the Hamiltonian [Z(J El| leading to a single 
band in the excitation spectrum. This is very similar to 
the feature displayed in Fig. [5] and Fig. [5] for small per- 
pendicular momenta k±. Related results have also been 
obtained in the Gross-Neveu model [43| and, recently, 
selfconsistent solutions have also been found for complex 
gap functions (li ]. 

The case of a superconductor with two non-degenerate 
fcrmion species in 1 + 1-dimensions has been discussed in 
Ref. It was found that there is a transition at 5 fx = 
^Abcs from the BCS phase to the solitonic phase, which 
persists to any larger value of Sfi at zero temperature. A 
selfconsistent solution of the gap function in this case is 
given by 

A 1+1 (x) = n\fv sn(/c(a; — x ); v) , (81) 

where sn(£ ; v) is a Jacobi elliptic function with elliptic 
modulus v 1 . The Jacobi elliptic function has the proper- 
ties sn(0; v) = 0, J|sn(£ ; ^)|g=o = 1) an d for < v < 1 
it is periodic in £ with period <iK(i>), where K(v) is the 
complete elliptic integral of the first kind. For the limiting 
cases, we have sn(£; 0) = sin£ and sn(£; 1) = tanh£, i.e., 
for v — > the gap function becomes sinusoidal, whereas 
in the limit v — > 1 we recover a single soliton. 

With Eq. (f8"Tj) . the ground state is obtained by min- 
imizing the thermodynamic potential in k and v. It is 
remarkable that - apart from xq, which is just an arbi- 
trary shift - the gap functions can be characterized by 
only two parameters, k and u, which can be related to 
the period AK(v)/n and to the slope n 2 ^/v at the zero 



crossings. In particular, the amplitude n^v is completely 
determined by these two quantities. Also note that k is 
just the slope at the zero crossings divided by the ampli- 
tude and we can therefore identify 2/k with the soliton 
size. Moreover, for constant K, if v is varied within an 
interval close to 1, the period changes strongly, whereas 
the amplitude and the slope stay nearly constant. This 
is very similar to our observation in Fig. [2] that the am- 
plitude and the shape of the transition region is almost 
independent of the period. 

One may thus wonder, whether the onc-dimensional 
gap functions in the 3+1 dimensional system are given 
by similar functions. Since, at least in the weakly coupled 
regime, all dynamics is constrained to a close shell around 
the Fermi surface, we could study the system on small 
patches at the Fermi surface, where the 3 -I- 1-dimensional 
problem can be reduced to 1 + 1-dimensional ones by ap- 
plying a quasi-classical approximation [49[ . The essential 
difference to the real 1 + 1 dimensional case is the fact 
that in 3+1 dimensions the momenta of the paired quarks 
are in general not parallel to the wave vector q of the con- 
densate. Therefore the quasi 1 + 1-dimensional problem 
needs to be solved for all directions [Hj]. However, 
projecting the wave-vector of the condensate on a given 
direction renders the relation between shape and period 
of the order parameter and its amplitude to depend on 
the direction. For this reason, the 1 + 1-dimensional ap- 
proach for obtaining a selfconsistent solution cannot be 
extended trivially to 3 + 1-dimcnsions. 

On the other hand, we already observed in the dis- 
cussion of Fig. [5] that, just as in 1 + 1 dimensions, the 
gap functions seem to depend on two independent scales, 
\q\ and A^cSj and that the latter determines both, the 
amplitude and the soliton size. Taking into account the 
complications discussed in the previous paragraph, this 
could be understood if the pairing is dominated by re- 
gions of the Fermi surface with a fixed azimuthal angle 
with respect to q. This idea is also inspired by Ginzburg- 
Landau investigations where the pairing mechanism of 
the fermions shows up more explicitly [34| . In particular 
for the FF phase at weak coupling, the pairing is concen- 
trated around rings on the Fermi surface where the rela- 
tive angle between k and g* is given by cos Off ~ j^- With 
this picture in mind we could imagine a similar regime 
of the Fermi ball to dominate the pairing in the general 
inhomogeneous phase with the order parameter varying 
in only one dimension. We would then expect that the 
amplitude of the gap function is again proportional to 
the inverse size of the soliton, but with a different pro- 
portionality constant than in the 1 + 1-dimensional case. 
Based on these arguments we suggest the fit 

A fit (z) = v4sn( K (z-z );^) (82) 

as a parameterization of our numerical results. To de- 
scribe a gap function with period and to comply with 
our convention of taking A(z) even with a maximum at 
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z — 0, we choose 
A = A(0) , 



4K(v)\q\ 



-'o 



7T 

"m 



(83) 



Then the only fit parameter left is v. As a measure for 
the quality of the fit we consider the relative deviation 



A - A 



fit 2 



A - A 



fit 2 



IA 



fit || 2 



A(0) 



vK{y) 



K(v)-E(v) 



,(84) 



in the parameterization. We have chosen Sfi = 0.7Abcs- 
However, we remind that the gap function is very insen- 
sitive under variation of 5[A, as pointed out in the context 
of Fig. [TTJ] It turns out that gap function is described re- 
markably well by the parameterization in Eq. (|82[) with 
only one fit parameter. The relative deviation in the L 2 - 
norm is of the order of one percent 7 . To illustrate this 
further, we present the "worst" case at |<fj = 0.6Abcs 
in Fig. [T2j The difference between the two functions is 
barely visible from plot. 



where E(y) is the complete elliptic integral of the second 
kind and II. II 2 is the L 2 -norm. 
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FIG. 13: cos 8 as defined in Eq. 
Sfi = 0.7A B cs- 
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FIG. 11: Relative error between the numerically determined 
gap function and its best fit by Eq. (|82|l as function of the 
period. Here we have chosen S/j, — 0.7Abcs- 
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FIG. 12: Numerical result for gap function A(z) compared 
to the best fit Afl t (z) at Sfi — 0.7Abcs and |<fj = O.GAbcs- 
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FIG. 14: Inverse soliton size k as function of |<fj for optimum 
fit to numerical result. 8fi = 0.7Abcs- 



Having found the astonishing agreement between the 
analytical form inspired from 1 + 1-dimcnsional investi- 



In Fig.[TT]we present the relative error of our numeri- 
cally obtained solutions compared to the best choice for v 



7 In view of the involved numerics and the marginal deviation we 
even do not want to exclude that the functions are identical. 
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gations and our numerical results, we are naturally lead 
to the assumption that the underlying dominant pair- 
ing mechanism is indeed dominated by fermions on the 
Fermi surface whose momenta have a fixed azimuthal an- 
gle measured from the direction of q. If this is true, we 
would expect that the amplitude of the gap function is 
again proportional to k, but with a different proportion- 
ality constant as in Eq. (f8"Tj) . 

To establish this connection, which has so far been ne- 
glected in our parameterization, we assume that the pair- 
ing in directions n with a particular azimuthal angle 9 to 
the direction of q is dominant. Furthermore, we assume 
that the gap function in 3 + 1 dimensions is just given by 
the 1 + 1 dimensional result, Eq. (f8"Tj) . if we go along a 
line in such a direction, i.e., 



where xg(z) 
given z and 



A 3+ iO) = A 1+1 {xg(z)), (85) 
is the coordinate in n direction for a 



cos 



Ai + i(x e ) = ngy/v sn(n e (xo - xefl)\v) (86) 

corresponds to Eq. (|8"T|). Hence, if we take k$ = kcos8, 
we obtain 



A 3 +i(z) 



\9\/v sn(ft(z — zq); iA . (87) 



Comparing this with our ansatz Eq. 
we get 



and Eq. 



cos 9 



A 



ttA(O) 



A^K( V )\q\ 



This quantity is shown in Fig. 1131 It turns out that our 
suggested picture breaks down at |<f| > 0.9Abcs, where 
cos 9 becomes larger than unity. This is related to the fact 
that the maximum of the gap function in our numerical 
results remains large also when its shape becomes more 
and more sinusoidal towards larger values for \q\. In the 
limiting case we have v — ► and Eq. (|88|) blows up. 

On the other hand, we find that cos# is almost con- 
stant for \q\ < 0.6 Abcs at a value of cos 6* w 0.62. This 
value corresponds to an opening angle 29 w 100°, some- 
what larger than the typical opening angle of 67° in the 
FF phase. Note, however, that the situation in the FF 
phase is rather different. There the pairing can be un- 
derstood most intuitively by shifting two Fermi spheres 
with radii \x + S/i and p, — Sfi by q in opposite directions 
[32j. In weak coupling, the physically relevant regime is 
then constrained around the intersection of the two Fermi 
surfaces, and one finds that cos 9 pp w rjjj-, as mentioned 
before. This means that 9pp is not constant as a function 
of |g| and only possible if \q\ > S/i. However, since in the 
FF phase \q\ by itself is almost constant ~ 0.9Abcs, one 
typically finds cos9pp 0.83. 

For the general one-dimensional inhomogeneous phase, 
this is quite different. In this case, as we have seen, \q\ 
varies strongly and can become arbitrarily small. Never- 
theless, according to Fig. [T31 the associated angle is ap- 
proximately constant, at least for \q\ < 0.6Abcs- This 



regime, where FF pairing would not be possible at all, 
agrees roughly with the realm of the soliton lattice, see 
Fig. [2 Here the gap functions are no longer dominated 
by the lowest Fourier components A^±i, see Eq. (|72|). 
but higher harmonics are important as well. In fact, as 
we have seen earlier, the scale \q\ which is related to the 
inverse distance of the solitons becomes rather irrelevant 
for the dynamics, which is more closely related to the 
inverse soliton size k. 

The latter is displayed in Fig. [TJ] It is interesting to 
see that its value is comparable with the corresponding 
scale in the FF and indeed almost independent of the 
periodicity of the solitons. From the numbers which can 
be read off from the figure we find that the soliton size 



^ is about 0.3 - 0.4 in units of 
agreement with Fig. [5] 



: . This is in good 



IV. SUMMARY AND CONCLUSIONS 

We have studied inhomogeneous pairing in relativistic 
imbalanced Fermi systems within a model of the Nambu - 
Jona-Lasinio type. The analysis was performed within 
mean-field approximation, but without restriction to the 
Ginzburg-Landau approach. In the set-up of the model 
we mainly focused on color superconducting phases, but 
the method is rather general and could be applied, e.g., 
to condensed matter physics or cold atomic gases as well. 

In Sec. [Til we developed the general formalism for cal- 
culating the thermodynamic potential. We considered a 
certain class of gap functions, which are time independent 
and periodic in space. The color and flavor structure of 
the gap functions was chosen according to the most im- 
portant pairing patterns in dense QCD, CFL and 2SC, 
but the extention to other phases is straight forward. 

We pointed out that the regularization of the theory 
needs to be addressed with special care in order to get 
a consistent description of homogeneous and inhomoge- 
neous phases. For instance, a straight-forward generaliza- 
tion of the three-momentum cutoff scheme to inhomoge- 
neous systems leads to undesired artifacts. To avoid these 
problems we suggest a Pauli-Villars-like regularization 
scheme, which can be derived via proper-time regular- 
ization. In this scheme the divergencies are regulated by 
restricting the free energies of the quasiparticles, rather 
than their momenta. As a consequence, the results are 
rather insensitive to the choice of the cutoff parameter if 
the gap in the homogeneous phase is fixed. In particular, 
model independent results in the weak-coupling limit are 
reproduced correctly. 

In Sec. Mil we discuss first numerical results obtained 
within this framework. To this end, we considered a sim- 
plified model with two fermion species at a chemical po- 
tential difference S/i. Basically, this corresponds to a re- 
striction to the 2SC phase in high-density approxima- 
tion. Moreover, in order to keep the numerical effort at a 
tractable level, we restricted ourselves to real gap func- 
tions with general periodic structures in one dimension. 
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With this ansatz we found that the inhomogeneous solu- 
tions are favored against homogeneous superconducting 
(BCS) and normal conducting phases in a Sfi-window 
which is about twice as wide as for a simple plane-wave 
ansatz (FF phase). The main effect is seen at the lower 
end of this window, i.e., towards the boundary to the 
BCS phase. In this region we observe the formation of a 
soliton lattice. With lowering the distance between the 
solitons increases and eventually diverges at the critical 
point. In this way the inhomogeneous phase is continu- 
ously connected to the BCS phase. On the other hand, at 
the upper end of the window where the inhomogeneous 
phase is favored, the gap functions are sinusoidal, and 
the transition to the normal conducting phase is of first 
order. 

It is worth noting that similar investigations of inho- 
mogeneous phases of cold atoms in the unitary regime 
using a local density functional approach (El, l53t | give ex- 
actly the opposite ordering of the phase transitions, i.e., 
first order for the transition from the homogeneous to 
the inhomogeneous phase and second order from the in- 
homogeneous to the normal phase. In our approach the 
behavior is however not unexpected: On the one hand 
side, mean-field investigations in 1 + 1-dimensions find 
the same behavior [4|| for the transition from the homo- 
geneous to the inhomogeneous phase. On the other hand, 
the phase transition near the tricritical point at finite 
temperature can be analyzed in a generalized Ginzburg- 
Landau approach [54], HH . The obtained results are per- 
fectly consistent with ours. At this point we would also 
like to mention that within the generalized Ginzburg- 
Landau approach the energetically preferred inhomoge- 
neous phase near the tricritical point has an order pa- 
rameter only varying in one dimension, as in our inves- 
tigation. It is however not clear whether this persists to 
zero temperature. 

We also studied the quasiparticle excitations in the in- 
homogeneous ground state as functions of the momentum 
k in the B.Z. As the inhomogeneous gap functions break 
the rotational symmetry, the spectra depend on both, the 
modulus and the direction of k, leading to an interesting 
band structure. We found that the spectrum is "almost 
gapped" with a few low-lying modes which correspond 
to the solitons. As a consequence, the gap functions are 
almost insensitive to variations of <5/z if the period is kept 
constant. 

Finally, we compared our solutions for the one- 
dimensional gap functions with the known analytical so- 
lutions of the corresponding 1+1 dimensional theory, i.e., 
Jacobi elliptic functions. We found that we can achieve 
an excellent fit in 3+1 dimensions, if we allow for one 
additional parameter. This can be interpreted as an ef- 
fective 1+1 dimensional behavior which comes about if 
the pairing is dominated by fermions on the Fermi sur- 
face with momenta at a fixed azimuthal angle relative to 

The present paper is only a first step towards a more 
complete description of inhomogeneous pairing in rela- 



tivistic systems. In the future, these studies should be 
extended in several directions: 

We have calculated the gap functions as functions of 
a given chemical potential difference <5/j. It would be in- 
teresting to see how this translates into density profiles 
of the different fermion species. As argued at the end of 
section MI D[ we expect that the density difference will 
be close to zero at the constant plateaus of the gap func- 
tions and that it will be peaked at the zero crossings, 
i.e., in the solitons. Having worked this out, one could 
apply this to describe more physical situations, like glob- 
ally neutral two-flavor quark matter in beta equilibrium 
or imbalanced atomic systems with fixed concentrations. 

The analysis could be extended to finite temperature. 
In the context of color superconductivity one should also 
study other pairing patterns, like the CFL phase. Even- 
tually, all this should be combined to obtain a phase dia- 
gram where homogeneous and inhomogeneous phases are 
treated on an equal footing. It would be interesting to see 
whether this can cure the problem of chromomagnetic in- 
stabilities in the phase diagram of neutral quark matter. 

Of course, it would be desirable to relax the restric- 
tion to one-dimensional gap functions and to study two 
and three-dimensional crystalline structures. In fact, ear- 
lier investigations of crystalline color superconductivity 
revealed that three-dimensional crystals are the most fa- 
vored inhomogeneous solutions [H, [H, H3| • These anal- 
yses, however, were performed in Ginzburg-Landau ap- 
proximation, which turned out to be particularly prob- 
lematic in the two-flavor case [34j. Moreover, the crystal 
structures were restricted to superpositions of a finite 
number of plane waves whose wave vectors all have the 
same length. This means that, e.g., the non-sinusoidal 
one-dimensional solutions which we found close to the 
BCS regime were not included. It is therefore not clear 
a priori which solutions are favored in a more com- 
plete analysis with arbitrary one-, two-, and three- 
dimensional periodic structures, in particular also be- 
cause a one-dimensional solution is expected near the tri- 
critical point [54j, [55j. In principle, this could be studied 
within our framework. In practice, of course, the numer- 
ical solutions will become very involved, and one has to 
see how far one can get. 
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APPENDIX A: GINZBURG-LANDAU 
EXPANSION AND REQUIREMENTS ON 
REGULARIZATION PROCEDURE 

Starting from the mean-field thermodynamic potential, 
Eq. (|20p , we can derive a Ginzburg-Landau functional by 
expanding the potential in powers of the pairing gap. To 
that end, we split the inverse propagator, Eq. (|19p . into 
a free part and a gap dependent part, 



S- 1 = Sq 1 + A . 



(Al) 



where So — S\ A _ and expand the logarithm in Eq. |2T 
in a power series in A: 



Trln 



Trln 



T 



A 



1 \ °° 1 



Since Sq A has only off-diagonal terms in Nambu-Gor'kov 
space, only even n contribute to the trace. To lowest non- 
trivial order, we thus get 







MF 



MF\ A=0 



A q k 



■0(A 4 ) 



(A3) 



For the evaluation of the trace, we need the explicit ex- 
pressions for 5*o and A. From Eq. ([TO)) we get 



So 



S+ 
5V 



and A 



A H 



(A4) 



with 



A 



Ik 



(A5) 



and A = — (A + )L Furthermore, we get from Eq. @ 

A qh =Y J A A , qk r A \ A , (A6) 

A 

while 5*0 is diagonal in flavor and color, 

5 ± =diag /c ((5 ± ) /c ) , (A7) 

with flavor indices / € {u, d, s} and color indices c £ 
{r, g, b}. Inserting these expressions into Eq. (|A3|) , the 
trace is readily evaluated. As a consequence of momen- 
tum conservation, we find that unequal Fourier compo- 
nents of the gap function do not interfere at this order. 
Likewise, there is no interference between two compo- 
nents with different color-flavor indices. We thus obtain 



an incoherent sum over \A A<gk | 2 , which can be combined 
with the sum in Eq. (|A3[) to get 



Qmf — ^mf|a=o 

4EE<i a ^i 2 + °( a4 )- ( A8 ) 

A q k 

The coefficients &A,q k are given by 



T 



(A2) , J_ 



2V ^ ^ [( S o)vr(Pn + Qk)l5 (5*0 )dg(Pn)l5 

+ tr D [(S + ) U g(p n + qk)lb {Sq )dr{Pn)l5 
+ tr D [{S£)dr(Pn + (lk)lb {So)ug(Pn)j5 
+ tr D [{S£) dg (pn +Qk)j5 (Sa)ur(p n ) T5 

1 

AH ' 



(A9) 



and analogously for a5. qk and o.7, qu - In the 2SC phase 
only A2 i9fc is nonvanishing, and we may drop the index A. 
The coefficients then essentially depend on the chemical 
potential difference <5/i, defined in Eq. (Ti^|) . 

The remaining traces in Dirac space, denoted by trjj, 
are trivial. As before, the sum over p n should be read as 
a Matsubara sum and a sum over the three-momentum 
p n , cf. Eq. HH). In the infinite volume limit, the latter 
should be replaced by an integral. Moreover, we have 
to introduce a regularization procedure to render this 
integral finite. 

As pointed out above, the quadratic term of the 
Ginzburg-Landau functional is an incoherent superposi- 
tion of the contributions from the different Fourier com- 
ponents. In particular, the coefficients a qk are completely 
independent of each other and should therefore be equal 
to the coefficients in the Fulde-Ferrell phase at given q 
and S/i. 

The determination of a{qk,Sfi) in the Fulde-Ferrell 
case is a well-known exercise (see e.g. Refs. [HI, |32j])- In 
a weak-coupling expansion at high densities one finds 



Ho) = -2(a{\q\,SiJ,) + 5a re g) , 



(A10) 



where 



a{q,5n) 



-1 



^ln 

2q 



5/j, 



A 2 

^BCS 



,9 - 

(All) 

while Sa reg depends explicitly on the regularization in 
the sense that the cutoff dependence of this term cannot 
be absorbed in an observable, like Abcs- As we pointed 
out in section Ul CI such terms should be avoided in order 
to keep the results free from regularization artifacts. In 
fact, in the regularization schemes usually employed for 
the Fulde-Ferrell phase [H], [HJ , the dependence on the 
regularization can completely be absorbed into a Abcs 
dependence, and Sa reg vanishes. 
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However, it is not obvious whether this regularization 
can be generalized to arbitrary gap functions in the ther- 
modynamic potential. Using a sharp cutoff A for in- and 
out-going momenta instead, we obtain 



Sa cutoff = 

reg 4(A 



klA 2 



m)m 2 



(A12) 



This excludes this most naive regularization scheme. For 
the proper-time regularization introduced in section lll CI 
on the other hand, we find 



r prope 
ulx reg 



= 



(A13) 



in weak coupling. We therefore suggest to use this scheme 
for a consistent regularization of inhomogeneous phases. 

Finally, we would like to add two comments: First, in- 
stead of starting from Eq. (|2Tjl . the Ginzburg-Landau 
functional could be derived equally well by expanding 
Eq. f25]). At T = this means that the Ginzburg- 
Landau expansion corresponds to a perturbative expan- 
sion of the eigenvalues of the Hamiltonian, with the pair- 
ing gap treated as perturbation. However, since some of 
the eigenvalues of the "unperturbed" Hamiltonian vanish 
at the Fermi surfaces of the particles, whereas the pair- 
ing gap is large, the perturbative expansion does actually 
not converge for the most interesting eigenvalues near the 
Fermi surfaces. As an example consider a BCS-like dis- 
persion relation, 



E(p) = V(P-M) 2 + |A| 2 
I I |A|2 



2|p- M | 



0(|A| 4 ). 



(A14) 



which obviously does not converge near p = fi. The 
Ginzburg-Landau expansion may therefore suffer a simi- 
lar problem. 

Second, the Ginzburg-Landau expansion may be useful 
to identify the divergencies in the thermodynamic poten- 
tial. In 3 + 1 dimensions it can be used to show that all 
divergencies can be absorbed by adding the expressions 
c 2 J |A(x)| 2 and C4 J \A(x)\ 4 with appropriate coefficients 
C2 and C4 to the thermodynamic potential. 



APPENDIX B: SYMMETRY OF THE 
EIGENVALUE SPECTRUM UNDER S^i 

In this appendix we discuss properties of the eigen- 
spectrum of the Hamiltonian Ha,s^ in Eq. (|5Tj) . For the 
case of a real gap function A (2) we would like to show, 
that in terms of the eigenvalue spectrum {Ex} of Ha,^ 
the eigenvalue spectrum of Ha.-s^ is given by {— Ex}. 
As a consequence, all eigenvalues of the total Hamiltonian 
Hhde, Eq. ([SO]), come in pairs {Ex, -Ex}. 

Restricting to real condensates, we can write 



Ha : *> = H a 3 + Acti - Sfj,l , 



(Bl) 



where (H )p m ,p n = (p m - l*)6p m ,f n , ^ m ,v~n = A Pm - P; 
and {<Ji} are the conventional Pauli matrices. Here we 
used that the Fourier components obey A qk — A*_ qk for 
a real gap function, see Eq. (TTj) . From {tr a , at,} = 2<5 a fcl 
we then obtain with J = ia% 



Ha.qJ = —JUa.o- 



(B2) 



Consequently an eigenvector v with eigenvalue E will al- 
ways come along with an eigenvector Jv with eigenvalue 
—E, i.e., eigenvalues of Ha,q come in pairs (Ex, —Ex). 

More specifically we can choose the eigensystem of J 
as a basis in which 



Ha,o = U 





-Ho - iA 



-Ho + iA 



with 



U 



_L (- 1 1 
V2 U 1 



Therefore we have 
with 



H I + V+ 



Hi + y_ 



V± = A 2 ±i{H ,A] 



, (B3) 



(B4) 



C/t , (B5) 



(B6) 



One can easily show that the spectra of the two operators 
Hq + V± are identical and given by {Ef}. This fact is 
sometimes referred to as being isospectral. The operators 
V± are furthermore very similar to potentials of fermions 
and bosons in supersymmetric quantum mechanics with 
superpotential A. 

Turning on again the chemical potential shift, we note 
that 



Ha. 



Ha,o — Sfil 



(B7) 



i.e., Sfj, just leads to a shift of the eigenvalues of Ha,o- 
Hence, for each pair (Ex,— Ex) of the eigenvalues of Ha,o, 
there is a corresponding pair (£\=F<5/i, —ExT^n) °f eigen- 
values of Ha,±&ii- Therefore we get the claimed connec- 
tion between the spectrum of Ha,8^ and Ha,-8^- 

Note that, in this proof, we made explicitly use of the 
fact that the matrix A is hermitian, which, in turn, was 
a consequence of our assumption that the gap function 
A(x) is real. This is the case for the inhomogeneous solu- 
tions discussed in Sec. IIII Cl and thereafter, as well as for 
the BCS phase. It is not true for the FF phase, where the 
gap function is complex. However, in this case we can ex- 
plicitly see from Eq. (|7Tj) that the eigenvalues E±(k) go 
over into —E^z(-k) when we replace 5/j, by — 6fx. There- 
fore, since we integrate over k in Eq. (|70[) . Ha,s^ and 
Ha,-s^i contribute equally to the thermodynamic poten- 
tial. 
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